Acorn User: China

home *** CD-ROM | disk | FTP | other *** search

/ Acorn User: China / Acorn User China CD-ROM (UK) (Disc B) / Acorn User China CD-ROM (UK) (Disc B).bin / STUTTGART / DATABASE / ULTIMATE / a252div.txt < prev next >

Wrap

Internet Message Format | 1993-08-09 | 20.0 KB

From: chughes@maths.tcd.ie (Conrad Hughes) Subject: Re: fast divide Date: 31 Oct 91 13:15:24 GMT DEFFNdiv(A,B,C,R):LOCALA%:[OPTpass:MOVS C,A,LSR#31:RSBNE A,A,#0 :ADDS C,C,A,LSL#16:MOV A,A,LSR#16:ADDS A,B,A:SUBCC A,A,B:] FORA%=1TO31:[OPTpass:ADCS C,C,C:ADCS A,B,A,LSL#1:SUBCC A,A,B:]NEXT [OPTpass:ADCS R,C,C:RSBCC R,R,#0:]=0 This is the BASIC assembler code to do the division.. Use it this way: FNdiv(reg_top,reg_bot,reg_misc,reg_res) while in the BASIC assembler.. It inlines the code. Register reg_res contains 32768*(reg_top/reg_bot); reg_bot _must_ be negative (if you don't know what sign it'll be, tack the following in at the beginning: CMP B,#0:RSBPL B,B,#0:RSBPL A,A,#0 ...) , and reg_top can be either sign.. If you need slightly higher accuracy (I think I can increase the accuracy by half a bit) and don't worry about registers, tell me. reg_res can be any register, including reg_top, reg_bot or reg_misc. If you want to turn it into a loop, do so.. But remember that the fundamental unit of of the routine takes 3 clock counts, a branch-loop takes five! Also, it'll be straightforward enough to turn into a subroutine.. If you need any clarification of it, ask me... Someone posted another version which is faster for a logarithmic distribution of numbers, with 0 and 1 occuring most frequently; I haven't got it to hand, but could dig it up if nobody else posts it. Conrad From: torbenm@diku.dk (Torben AEgidius Mogensen) Date: 31 Oct 91 15:34:21 GMT Subject: Re: fast divide Sender: torbenm@freke.diku.dk Here is one I cooked up. It is from memory, so there might be slight errors. The comparisons before the first SUB ensure that no overflow happens on the ASL. Worst case time is 5 cycles per bit. Arguments in p,q. on exit, r = p div q, p = p mod q .div MOV r,#0 CMP p,q BLT nodiv CMP p,q ASL #1 BLT div0 CMP p,q ASL #2 BLT div1 ... \ repeat for 3..30 CMP p,q ASL #31 SUBGE p,q ASL #31 ADDGE r,#2^31 CMP p,q ASL #30 .div30 SUBGE p,q ASL #30 ADDGE r,#2^30 CMP p,q ASL #29 .div29 ... \ repeat for 29..1 .div0 SUBGE p,q ASL #0 \ equiv. to SUBGE p,q ADDGE r,#2^0 .nodiv MOV PC,Link From: gcwillia@daisy.waterloo.edu (Graeme Williams) Date: 2 Nov 91 20:13:28 GMT Subject: Re: fast divide Organization: University of Waterloo Ok, here is the fastest unsigned divide routine I know of. The 3 instruction loop (unrolled) wasn't my idea - but the piece of code that figures out whereabouts to jump into the unrolled loop is. On entry d and n should be both +ve. CMP d,#0 : BEQ end% ; Remove zero divisors CMP d,n : MOVGT m,n : MOVGT n,#0 : BGT end% ; Exit if denominator > numerator ; This next bit of code figures out roughly how many bits the quotient is ; going to have and hence how many times you need to run through the ; divide loop proper - its only worth finding out to the nearest 4 times .start% MOV cnt,#28 : MOV m,n,LSR#4 CMP d,m,LSR#12 : SUBLE cnt,cnt,#16 : MOVLE m,m,LSR#16 CMP d,m,LSR#4 : SUBLE cnt,cnt,#8 : MOVLE m,m,LSR#8 CMP d,m : SUBLE cnt,cnt,#4 : MOVLE m,m,LSR#4 MOV n,n,LSL cnt ADDS n,n,n : RSB d,d,#0 ; Have to flip the sign of d for the divide loop ADD cnt,cnt,cnt,LSL#1 ; Multiply cnt by 3 ADD pc,pc,cnt,LSL#2 ; Jump cnt instructions MOV R0,R0 ; Dummy instruction to take care of pipelining 32 copies of { ADCS m,d,m,LSL#1 : SUBCC m,m,d : ADCS n,n,n } .end% The best times occur when the numerator and denominator of your division are of similar size. The worst times when the numerator is very large and the denominator small. (The case where the denominator is larger than the numerator is killed off very early.) Suppose one has a distribution of numbers wherein the highest bit of the numbers is equally distributed between bit 0 and bit 31, that is, there is the same number of numbers between 2^n and 2^(n+1) for every n. (I think this is actually the case for most purposes) Now if a numerator and denominator are selected at random from this distribution subject only to the requirement that denom <= numer then the above routine will execute on average in 58.25 cycles (from the point start% in the code). The worst case is 116 cycles (top bit of quotient is bit 28 or higher) the best case is 32 cycles (top bit of quotient is bit 0 to 3) The inbetween cases get 12 cycles worse for every 4 extra bits in the quotient. Note that the average case isn't the simple average of the various cases but has to be weighted toward the better cases - this is because there are more ways (27) for numer and denom to differ by say 5 bits than there are (7) for them to differ by say 25 bits. Have fun. Graeme Williams From: zrzm0370@rusmv1.rus.uni-stuttgart.de (Joerg Scheurich) Subject: integer division Date: 10 Jun 91 11:25:36 GMT Hi So i know, there is no Assembler-command for integer-division in the ARM. Therefor a fast division-routine is a thing of common interest. In a earlier Posting i wrote about the bad quick hacked division-routine for the demo of my formula-tool (math. formula ==> assemblermacros). Now i wrote a better routine, but i don't know, if there is a better algorithm. The shorter and slower 32-Bit/32-Bit-division-routine requires 382 Cycles in 100 Bytes, with usage of 4 32-Bit-registers and 4 Byte memory. With unrolling the loop, the faster and longer 32-Bit/32-Bit-division-routine requires 182 Cycles in 692 Bytes, with usage of 4 32-Bit-registers. This is very faster then the old routine (which requires 2**33 Cycles in worst case and requires all time of the universe, if there is a division by zero ... ) Cause the INTERNAL-assembler-command of the INTEL 8086 for a 32-Bit/16-Bit-division requires 165-184 Cycles in 2-4 Bytes, with usage of 3 16-Bit-registers, i think, the result is ok. The result is nothing against the INTERNAL-assembler-command of the ARM for a 32-Bit*32-Bit-multiplication which requires 17 Cycles in 4 Bytes, with usage of 1-3 32-Bit-registers, but my Professor in digital electronics used to say: "If a program contain much divisions, the programmer is a idiot" ... Know someone a better (faster or shorter) algorithm ? so long MUFTI ( internet: zrzm0111@helpdesk.rus.uni-stuttgart.de ( from janet: zrzm0111%de.uni-stuttgart.rus.helpdesk@NFS-RELAY ) bitnet: ZRZM AT DS0RUS1I ) div-fast: --------- ; division routine ; 1. calculate flag if result is negate and convert operands to positiv ; 2. second "divide" with unsigned suczessive Approximation ; 3. fix the sign of result .MACRO INTdivide LDR R2,%2 LDR R3,%3 MOV R0,#0 CMP R2,#0 RSBLT R2,R2,#0 SUBLT R0,R0,#1 CMP R3,#0 RSBLT R3,R3,#0 MVNLT R0,R0 MOV R1,#0 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l2 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l3 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l4 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l5 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l6 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l7 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l8 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l9 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l10 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l11 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l12 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l13 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l14 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l15 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l16 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l17 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l18 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l19 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l20 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l21 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l22 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l23 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l24 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l25 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l26 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l27 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l28 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l29 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l30 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l31 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 \l32 ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 CMP R0,#0 RSBNE R2,R2,#0 STR R2,%1 .ENDM div-short: ---------- ; division routine ; 1. calculate flag if result is negate and convert operands to positiv ; 2. second "divide" with unsigned suczessive Approximation ; 3. fix the sign of result .MACRO INTdivide LDR R2,%2 LDR R3,%3 MOV R0,#0 CMP R2,#0 RSBLT R2,R2,#0 SUBLT R0,R0,#1 CMP R3,#0 RSBLT R3,R3,#0 MVNLT R0,R0 STR R0,minusflag MOV R0,#32. MOV R1,#0 \loop ADDS R2,R2,R2 ADCS R1,R1,R1 CMP R1,R3 SUBGE R1,R1,R3 ADDGE R2,R2,#1 SUB R0,R0,#1 CMP R0,#0 BNE \loop LDR R0,minusflag CMP R0,#0 RSBNE R2,R2,#0 STR R2,%1 .ENDM From: chughes@maths.tcd.ie (Conrad Hughes) Subject: Re: integer division Date: 11 Jun 91 14:10:31 GMT zrzm0370@rusmv1.rus.uni-stuttgart.de (Joerg Scheurich) writes: >382 Cycles in 100 Bytes, with usage of 4 32-Bit-registers >and 4 Byte memory. >With unrolling the loop, >the faster and longer 32-Bit/32-Bit-division-routine requires >182 Cycles in 692 Bytes, with usage of 4 32-Bit-registers. This divides a 64bit number by a 32bit one, with 32bit result, in 97 instructions/388 bytes/97 clock cycles.. Entry: A=MSW of 64bit no., B=LSW of 64bit no., C=32bit no. C _must_ be negative (if it's positive, RSB C,C,#0), A _must_ be positive (if not, RSBS B,B,#0:RSC A,A,#0) ADDS B,B,B ADCS A,C,A,LSL#1:SUBCC A,A,C:ADCS B,B,B [repeat these three 32 times] Exit: A=Remainder of division, B=result, C unchanged. ..if you don't need the remainder, get rid of the last SUBCC A,A,C. The result will be unsigned, and flags will (obviously) reflect the last addition. Adjust as you will to get 32/32 division (Initialise A to 0, or AB = numerator shifted some way) or any other variety you like.. Rolling up the loop it could get slightly complicated since every instruction requires flag information from the previous one ;-) Something like SUB count,count,#1:TEQ count,#0:BNE loop with count initialised to #32 might do the trick, but it'll slow it down to at least 288 cycles (so what if it only uses 32 bytes?? ;-}), which is pretty unacceptable unless you're on an ARM3.. This is off the top of my head and I'm nowhere near the Arc, so if it doesn't work (it should - I've typed it in enough times now..) flame me.. Conrad From: gcwillia@daisy.waterloo.edu (Graeme Williams) Subject: Re: really FAST integer division Summary: Avg. exec. time < 60 cycles - the fastest div in the West? Date: 28 Jun 91 19:07:58 GMT Organization: University of Waterloo Following several posts regarding division routines, in particular one with a lovely 3 instr. central loop I modified my old division routine into the following (supersonic, lightning fast, ZR rated, :-) ) 32bit routine: ;Essentially what's happening here is we are figuring out how far left ;we can shift the numerator before any actual dividing begins and thus save ;a significant number of useless trips through the (unrolled) 3 instr. loop. ; - Note we only go to the nearest 4 bits, trying to save another 2 bits ; results in 3 extra instrs. for a saving of 6 instrs. half the time ; and 0 instrs. the other half - i.e. no nett benefit. MOV cnt,#28 : MOV mod,num,LSR#4 CMP den,mod,LSR#12 : SUBLE cnt,cnt,#16 : MOVLE mod,mod,LSR#16 CMP den,mod,LSR#4 : SUBLE cnt,cnt,#8 : MOVLE mod,mod,LSR#8 CMP den,mod : SUBLE cnt,cnt,#4 : MOVLE mod,mod,LSR#4 MOV num,num,LSL cnt : RSB den,den,#0 ADDS num,num,num ;Now skip over cnt copies of the 3 instr. loop. ADD cnt,cnt,cnt,LSL#1 : ADD pc,pc,cnt,LSL#2 : MOV R0,R0 {ADCS mod,den,mod,LSL#1 : SUBLO mod,mod,den : ADCS num,num,num} x 32 copies The routine behaves as follows: On Entry: num - is a signed numerator A (but +ve) den - is a signed denominator B (but non-zero and +ve) On Exit: mod - contains A MOD B num - contains A DIV B The routine requires 4 registers and 452 bytes of memory. Execution time: Given numbers A and B (both +ve) chosen at random from a distribution in which the most significant bit of the numbers is uniformly distributed between bits 0 and 31 (i.e. the density of numbers goes as roughly log n) then the average execution time is: For A < B : 32 cycles (no division actually required) For A >= B: 58.25 cycles The best and worst cases are 32 cycles and 116 cycles respectively. Have fun. (Stuff like this ought to be worth good money!) Graeme Williams gcwillia@daisy.waterloo.edu P.S. With any luck there's only a couple of typos :-). From: dseal@armltd.uucp (David Seal) Subject: Re: Dividing Date: 5 Mar 91 19:47:12 GMT In article <8828@castle.ed.ac.uk> ecwu61@castle.ed.ac.uk (R Renwick) writes: > Can anyone help me with this problem: >In my super-duper, wire-frame model animating 'C' program, I represent >floating point numbers as 4-byte integers. This I do by shifting the number >being represented by 15 bits (multiply by 32768) in order that I can keep a >track of the fractional part. The reason for doing this is that I want >to make sure my program runs as fast as possible by not using floating >point arithmetic... > > What I want to do is do a divide a number by another number >and not lose alot of the precision. > The way I do this at the moment is by > result=((a/b)<<15) >But if a and b are close together then (a/b) loses alot of precison >(since the result is an int). > >Does anyone know how I could do an integer division without losing >precision. I cannot use result=(((a<<8)/b)<<7) or any other tricks >because a may be close to 32768<<15 already... >Perhaps if I can find the reciprocal of b in the above representation >and multiply it by a then I'll not lose the as much precision, but I >don't know how to do this (find the reciprocal) :-( Two quick suggestions: (1) The basic problem is that integer division stops when it determines the bit in the units place, whereas you want it to go on and calculate another 15 fractional binary places. So rather than using C's integer divide operator, write a long division routine using subtractions and shifts and force it to go on 15 iterations longer than normal. (This won't be much slower than C's integer divide, which essentially uses a long division routine anyway). If speed is very important, you'll probably want to either (a) write a well-optimised assembler routine; or (b) write it as a C macro, so that the compiler includes it inline. In the latter case, you may want to look carefully at the assembler produced by the compiler for an instance of your macro and experiment a bit - careful choice of the details of the long division algorithm can make a significant difference to the speed of the division. Incidentally, you may want to watch out for significant bits being lost during the last 15 iterations: this is perfectly possible and means that your result has overflowed. (2) Use a Newton-Raphson iteration to find the reciprocal of your divisor and then multiply - this presupposes that you've got an efficient way to do multiplications, both for the Newton-Raphson iteration and for the final multiplication. The iteration to find the reciprocal of a number A is (in pseudo-code - you'll have to translate to C): newx = initial approximation to reciprocal; repeat x = newx; newx = x * (2 - A*x); until x and newx sufficiently close together; The main thing to beware of in this is that your initial approximation must be reasonably close to the reciprocal - it must be greater than 0 and less than twice the reciprocal. Such an approximation can be found by finding the most significant 1 in A and choosing a value based on this - e.g. if A lies in the range 1.0 <= A < 2.0, its most significant 1 is at bit position 16, and a suitable starting point for the iteration would be 0.75. This gives us: Most sign. bit position Range for A Suitable starting point ----------------------------------------------------------------------- : : 9 2^(-7) <= A < 2^(-6) 1.5*2^6 10 2^(-6) <= A < 2^(-5) 1.5*2^5 11 2^(-5) <= A < 2^(-4) 1.5*2^4 12 2^(-4) <= A < 2^(-3) 1.5*2^3 13 2^(-3) <= A < 2^(-2) 1.5*2^2 14 2^(-2) <= A < 2^(-1) 1.5*2^1 15 2^(-1) <= A < 2^0 1.5*2^0 16 2^0 <= A < 2^1 1.5*2^(-1) 17 2^1 <= A < 2^2 1.5*2^(-2) 18 2^2 <= A < 2^3 1.5*2^(-3) 19 2^3 <= A < 2^4 1.5*2^(-4) 20 2^4 <= A < 2^5 1.5*2^(-5) 21 2^5 <= A < 2^6 1.5*2^(-6) 22 2^6 <= A < 2^7 1.5*2^(-7) : : Watch out for very big values - e.g. 2^16 is representable in this number system, but 2^(-16) isn't! Hope these suggestions help. David Seal dseal@acorn.co.uk or dseal@armltd.uucp